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The problem of object restoration in the case of spatially incoherent illumination 
is considered. A regularized solution to the inverse problem is obtained through a 
probabilistic approach, and a numerical algorithm based on the statistical analysis of 
the noisy data is presented. Particular emphasis is placed on the question of the posi- 
tivity constraint, which is incorporated into the probabilistically regularized solution 
by means of a quadratic programming technique. Numerical examples illustrating 
the main steps of the algorithm are also given. 



The inverse problem in optics consists of recovering the object by starting from its image. 
It can be regarded as a backward channel communication problem: messages can be conveyed 
back from the data set (the image) to reconstruct the signal (the object). The length of these 
messages is limited by the noise affecting the imaging process. This fact can be viewed as the 
necessity of truncating the eigenfunction expansions associated with the Fredholm integral 
equation, which gives the mathematical formulation of the inverse problem in optics. The 
latter can be formulated as follows [1]: Consider a one-dimensional object, illuminated by 
spatially incoherent radiation and imaged by a perfect optical instrument (i.e., without focus 
error) with a rectangular aperture. If we use f(x) to denote the spatial radiance distribution 
in the object plane, then the noise free spatial radiance distribution in the image plane is 
given by 
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R being the Rayleigh resolution distance. Here we have assumed finite extent of the object 
and the linear magnification of the optical instrument to be +1. 

When g(y) is given in the geometrical region of the image plane, i.e., the interval [—1, 1], 
the problem of object restoration is equivalent to solving the Fredholm integral equation 
of the first kind Af = g, where A is a symmetric compact positive definite operator. The 
solution of the integral equation is unique, i.e., the equation Af = has only the solution 
/ = 0. 

A formal solution to the equation Af = g can be given by means of an eigenfunction 
expansion, i.e., 

oo 

/w=Ef^)' (3) 



where g^ = (g,ipk) [(•>•) denoting the scalar product in L 2 (— 1,1)], and ipk(x) and are, 
respectively, the eigenf unctions and the eigenvalues associated with the operator A. If we 
add to the image g a perturbation such as the one produced by the noise measurement, then 
the equation Af + n = g = g + n, n(y) being a function describing the noise or the ex- 
perimental error, generally has no solution. Accordingly, the expansion Y^^iidk/^kj^kix), 
[9k — (g, i>k)] diverges if g does not belong to the range of the operator A. We are faced 
with the pathology of the ill-posedness of the problem, in the sense of Hadamard Q]: small 
perturbations of the data produce wide oscillations of the solutions. The problem requires 
regularization. One of the most popular methods in use is due to Tikhonov 0, 0] and con- 
sists of restricting the solution space by imposing suitable global bounds on the solutions. 
Within the framework of the eigenfunction expansions, Tikhonov's regularization provides 
a criterion for a suitable truncation of the series. Furthermore, the function f(x), which 
represents a spatial radiance distribution, is a nonnegative function, and, consequently, a 
nonnegativity constraint must be necessarily added in the mathematical formulation of the 
problem. 

At this point it is worth noting that, although the L 2 space is the most natural ambient 
when the method of the eigenfunction expansion is used and the Tikhonov's regularization 
is adopted, the L 1 space could present several advantages, as we explain below. First, the L 1 
norm of the intensity is the energy radiated by the object, and therefore an equality of the 
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following type, II/Hl^-i.i) 

= f_ 1 \f(x)\dx = E, (E is the energy radiated by the object), can 
clearly be interpreted from a physical viewpoint. Analogously, the quantity J, \ Af — g\dy 
represents the statistical fluctuations of the energy of the image. Finally, the positivity 
constraint is sufficient to restore the continuity of the inverse problem in the topology induced 
by the L 1 norm (see Appendix A). Conversely, the positivity is not sufficient to regularize 
the problem in the L 2 topology. All these considerations point towards choosing the L 1 norm 
in the mathematical formulation of the problem. But it can easily be shown, by exhibiting 
explicit examples, that the continuity restored in the L 1 topology is far from implying a 
continuity in the topology induced by the sup norm (see Appendix A). Moreover, it should 
be observed that the positivity constraint is sufficient to restore the continuity in the L 1 norm 
but not in the L 2 norm, simply because the L 2 convergence is more robust [note that f(x) has 
compact support] and represents precisely the minimal requirement in the actual numerical 
calculations. In view of this last consideration we are then led to choose the L 2 norm in the 
mathematical formulation of the problem, and, consequently, the eigenfunction expansions 
can be used. However, we shall not regularize the problem by using Tikhonov's variational 
method, which requires a bound of the type: ||/||l 2 (-i,i) ^ M, (M is constant), whose 
interpretation is not clear from the physical viewpoint. Instead, we work out the problem 
with a probabilistic approach, which makes it possible to split the Fourier coefficients of 
the noisy data function [i.e., = (g,ipk)} into two classes: one comprises those coefficients 
from which a significant amount of information can be extracted; the other, those Fourier 
coefficients that can be regarded as random numbers because the noise prevails on the 
information content. We can thus construct an approximation that satisfies the requirements 
of the probabilistic regularization as explained in Subsection III Al Questions related to the 
actual numerical computation of the probabilistically regularized solution are discussed in 
Subsection III Bl where an algorithm based on the analysis of the autocorrelation function of 
the noisy data is presented. However, at this stage, the positivity constraint still remains to 
be satisfied. This point is addressed in Section IIII[ where the numerical issues concerning 
the construction of a nonnegative approximation are discussed, and examples of numerical 
tests are also given. Finally, in Appendix A, the ill-posedness of the problem in various 
topologies and the role played by the positivity constraint will be illustrated. 
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II. PROBABILISTIC REGULARIZATION AND STATISTICAL METHODS 

A. Probabilistic Regularization 

As remarked in SectionHJ the ill-posed character of the inverse problem is derived precisely 
from the fact that the data function g(y) is corrupted by the noise that, hereafter, will be 
represented by a bounded function n(y), which is supposed to be integrable in the interval 
[— 1, 1]; i.e., g = g + n. Furthermore, we set 

\\9-g\\m-iA) = IMIl2(-i,d ^ e - ( 4 ) 

In general, the perturbation produced by the noise is such that the noisy data function g 
does not belong to the range of the integral operator A, and, consequently, the eigenfunc- 
tion expansion YlkLiiSk/ '^k)i>k generally diverges. Further, even assuming that the noise 
perturbation is gentle enough for g to belong to the range of the operator A, in this case the 
noise n still precludes us from exactly knowing the noiseless data function g. If we have two 
distinct data functions g^ and g^ 2 \ we may not attribute to them two distinct solutions 
and /( 2 ), if the data distance is less than e, i.e., if \\g^ — g^ ||l 2 (-i,i) ^ e - Accordingly, 
we can consider to be distinguishable only those data such that \\g^ — g^ \\l 2 (-i,i) > £■ In 
conclusion, even if the series YlkLi{9k/ *k)ipk converges in the L 2 norm, it is meaningless to 
push the eigenfunction expansion beyond a certain value ko (i.e., a certain truncation point) 
that depends on e. 

The most intuitive and simple truncation method yielding a regularized solution consists 
of writing an approximation of the following type, 

fo(x) = V] T^fcO), (5) 

— Ah 

k=l K 

where k (e) is the largest integer such that ^ e [assuming, in this particular case, the a 
priori bound that the set of the input signals belongs to the unit ball in L 2 (— 1, 1)]. In fact, 
it can be proved □ DD that / (a;) converges weakly to f(x) even if g does not belong to 
the range of the operator A, i.e., 

hm(f-f ,v) L2{ _ ltl) =0 , [\/veL 2 (-l,l), IMI^-n) ^ 1]. (6) 

. Remark: If g G range (A), then it can be proved by use of the Kolmogorov e-entropy theory 
8| that fco(e) is strictly related to the maximal length of the messages L max (e) that can be 
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conveyed back through the communication channel associated to the operator A; in fact, we 
can prove that, for sufficiently small e (see also, Ref. 0), 

L m Ue)^2 ko{€)lo ^ 1/e) . (7) 
But the approximation fa(x) presents several defects: 

1. The solutions must be restricted to a bounded subset such as the unit ball in L 2 (— 1, 1) 
or, equivalently, to a bounded subset of the type: ||/||x^(-i i) ^ M, M = constant, 
whose physical interpretation is not transparent, as noted in Section |l] 

2. Only a weak convergence of fo(x) to f(x) is guaranteed, whereas at least the L 2 -norm 
convergence should be required for practical applications. In addition, fo(x) does not 
generally satisfy the positivity constraint. 

3. The truncation criterion ^ e (or ^ c/M) does not guarantee that the approxi- 
mation fo(x) really does pick out the Fourier components of the noisy data that are 
likely to carry exploitable information about the unknown solution and, at the same 
time, reject the ones dominated by the noise. 

To overcome all these difficulties, we turn the problem in a probabilistic form. With this 
in mind we rewrite integral equation JIJ in the following form: 

M + C = v, (8) 

respectively, are Gaussian weak random 



where £, ( and rj, which correspond to /, n and g 
variables in the Hilbert space L 2 (— 1, 1) (see Ref. llOr ) . A Gaussian weak random variable 
is uniquely defined by its mean element and its covariance operator; in the present case we 
use R,££, R^, and R m to denote the covariance operators of £, ( and r], respectively. 
Next, we make the following assumptions: 

I. £ and £ have zero mean; i.e. mg = = 0. 

II. £ and £ are uncorrelated, i.e. = 0. 

III. -R^ 1 exists. 
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The third assumption is the mathematical formulation of the fact that all the components of 



the data function are affected by noise. As shown by Franklin [see formula (3.11) of Ref. Ill| . 
if both signal and noise satisfy assumptions (I) and (II), then R m = AR^A*+R^, where A* 
indicates the adjoint of A, and the cross-covariance operator is given by: R^ v = R^A*. We 
also assume that R^ depends on a parameter e that tends to zero when the noise vanishes, 
i.e., R^ = e 2 N, where N is a given operator (e.g., N = I for white noise). 

At this point, we turn Eq. (jHJ) into an infinite sequence of one-dimensional equations by 
means of orthogonal projections 



h£k + (k = Vk, (fc = l,2, ...), 



(9) 



where = (£, "0ifc)> Ck = (CiV'fc)) Vk = (v^k) are Gaussian random variables. Next, we 
introduce the variances p\ = (R^ip k ,ip k ), e 2 z/| = (R (( ip k ,ip k ), \\p\ + t 2 v\ = (R m ip k ,ip k ), 
without assuming that the Fourier components of £ (and analogously also for ( k and rj k ) 
are mutually uncorrelated. In view of assumptions (I) and (III) the following probability 
densities for £ k and ( k can be assumed: 



2tt p k 
1 



27T€V k 
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(fc = l,2,...), 



(A; = 1,2,...). 



(10) 
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Using Eq. Q, we can also introduce the conditional probability density p Vk (y\x) of the 
random variable r] k for fixed £ k = x, which reads as 



Pvkivl 
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x 



2tt ev k 
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(12) 



'2ir ev k 

Let us now apply the Bayes formula that provides the conditional probability density of £ k 
given r) k through the following expression: 



Ptk( x \y) 



Pm(v) 



Thus, if a realization of the random variable r] k is given by g k , formula (JT3j) becomes 

2 



Pt k ( x \9k) = A k exp 



x 



2p 



exp 



2e 2 v 2 



x — 



9k 



(13) 



(14) 



Now, the amount of information on the variable which is contained in the variable r] k , 



can be evaluated. We have 



12 



J(£k,Vk) 



ln(l 



k)i 



(15) 



7 



where 



r 



k 





E{^77fc} 
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E{|6I 


2 }E{k| 2 } 



(16) 



Thus, we obtain 



'(&,w) = ih(l + ^|)- (17) 

Therefore, if A^p^ < ez/^, then J(£k,Vk) < (hi2)/2. Thus we are naturally led to introduce 
the following sets: 

I = {k : AfcPfc > ez/ fc } , J\[ = {k : \ k p k < eu k } . (18) 

If we revert to the conditional probability density (JHJ), this can be regarded as the 
product of two Gaussian probability densities: pi(x) = A k ^ exp (—x 2 /2pl) and p2(x) = 
A k 2 ^ exp {-(A 2 /2e 2 z/ 2 ) [x - (<7fc/A&)] 2 }, (A k = A k ■ A k 2 ^), whose variances are given by p 2 
and (ez/fc/A/c) 2 , respectively. Let us note that, if k 6 X, the variance associated with the 
density p2{x) is smaller than the corresponding variance of p\(x) and vice versa if k G A/". 
Therefore it is reasonable to consider as an acceptable approximation of the mean value 
given by the density p2(x) if k e X, but the mean value given by the density px(x) if k e Af. 
We can then write the following approximation: 

{q k /X k k E I 
9kl (19) 
keM 

Consequently, given the value g of the weak random variable T], we are led to the following 
estimate of £, which, using the notation of Ref. y, reads as 

Bg = iN*- ( 2 °) 

kex Xk 

Next, consider the global mean square error e|||£ — -Br]|| 2 | associated with the operator 
B, introduced in (|20|). We can now prove the following theorem: 

Theorem 1 If the covariance operator is of trace class, and furthermore 

lim/ c ^ 00 (Afcp/ c /z// c ) = 0, then the following limit holds true: 



lim E< 



'{ll^-^|| 2 }=0. (21) 
Proof: See Ref. Q. □ 
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B. Statistical Analysis of the Noisy Data 

The application of the results achieved in Subsection III Al calls for statistical tools able to 
determine the two sets X and N '. In this section this issue is discussed, and the basic steps 
of a numerical algorithm for constructing the regularized solution Bg from the noisy data g 
are outlined. 

Splitting the Fourier coefficients into the sets I and M can be performed by computing 
the correlation function of the random variables r] k , i.e., the probabilistic counterpart of the 
coefficients g k , 

A /. . N _ gjfal - E {^fcl})(^2 - E i)fe})} 

v[k1 ' 2) E{( Vkl -E{ Vkl m^ E{( Vk2 -E{ Vk2 }yy/f {ZZ) 

In practice, only a finite realization {g k }i of the random variables rj k is usually available, 
from which estimates Sg of the autocorrelations can be obtained by regarding the data 
{9k}\ as a finite length record of a wide-sense-stationary random normal series jisl ]. In 
principle, the assumption of stationarity of the series {f]k} is not strictly true, because the 
moments of the random variables i]k generally depend on k, but, from the practical point of 
view, this is usually the only possible option. However, the stationarity assumption can be 
removed whenever estimates of ensemble averages of the series {rj k } can be computed. Thus, 
by recalling that the are normally distributed, by introducing the working hypothesis 



that the process {rjk} is stationary in wide sense [l4|, i.e., A T7 (/ci,/c 2 ) = A v (ki — k 2 ), and 
by assuming that the ensemble contains no strictly stationary subensembles that occur 
with probability other than zero or one, we can compute estimates of the autocorrelation 
coefficients by means of the ergodic hypothesis Q| equating ensemble and time (i.e., the 
index k in our case) averages. 

Among the numerous estimators of the autocorrelation function , one which is widely 
used by statisticians is given by 

N—n 



^2(9k - (9k))(gk+n - (gk+n)) 

fc=l 

N-n 

^2(9k - (g~k)) 2 ^2(gk+ n - ((jk+n)) 



6- g (n) = ^ n = 0,...,A-l, (23) 

' N-n N-n 1 - 



k=l k=l 




where 



. N—n . N—n 

^) = T7— 7 ^ (9k+n) = W — 9k+n- (24) 



N-n^ ai WT/ N-n 

k=l k=l 
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Equation (|23|). which is based on the scatter diagram of g~k+n against g k for k = 1, .., N — n, 
represents the maximum-likelihood estimate of the autocorrelation coefficients of two ran- 
dom variables rjk and whose joint probability distribution function is bivariate normal. 

To identify the structure of the series {g~k}i so as to separate the correlated components 
from the the random ones, it is necessary to test whether 5g(n) is effectively zero. This 
question has been extensively discussed in Ref. f5j. Here, we briefly report on the main points. 
First we assume that there exists an index n such that for n > n , A v (ki — k 2 ) = A 7? (n) will 
vanish. This index n is actually recovered recursively as follows: the series {(jk}i is first 
supposed to be purely random, i.e., no = 0, the standard error a$(n;Q) is computed, and 
the smallest index n > such that |<5§(n)| > 1.96 a$(n; 0) is searched for. If such an index n 
is found, it becomes the new candidate n , i.e., we set n = n, and the whole procedure is 
repeated until no new index n is found. The large-lag standard error cr^n; n ) is evaluated 
by using the following formula, due to Bartlett [16]: 



1 



Formally, uq can be defined as 



N-n 



no 



v=l 



1/2 

, for n > n . (25) 



n = max{n ^ : Vn € (n, iV - 1] , | 6 3 (n) \< 1.96 a s (n,n)}. (26) 

Accordingly, the set Q of the lags corresponding to autocorrelation values effectively different 
from zero is defined as 

Q = {0 < n ^ n : | 5- g (n) |> 1.96 a s {n, 0)}. (27) 

Let Nq be the cardinality of Q. From Q we can construct Nq families Fi of pairs of Fourier 
coefficients defined by 

F l = {(g kl ,g kt+n M% nt \ niEQ, i = l,...,N Q , (28) 

from which the couples of coefficients that are likely to be correlated should be selected. 
At this point, the Fourier coefficients that are correlated are determined in a unique way by 
means of the following heuristic criterion: for any rii G Q, i = 1, ...,Nq, we select the pair 
(gk*,gk*+m) giving the maximum contribution to the corresponding autocorrelation estimate 
Sg(rii), i.e., we define fc* as 

fc* = arg max {\g k g k + ni \}, i = l,...,N Q . (29) 

£:6 1,7V— nJ 
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Accordingly, we can define the set of frequencies X that exhibit correlated Fourier coefficients 

I = {k$}pU + (30) 
where each element of X is counted only once. Finally, we can construct the approximation 

h(x) = J2yMx). (31) 

In theory, that is, for N —>■ oo, the Nj elements ki G X and the Nq numbers n ; £ Q 
are mutually constrained. In fact, any two coefficients k a ,k/3 G X must satisfy the pair- 
wise compatibility constraint requiring \k a — kp\ G Q. Moreover, it easy to see that 
the number N% of the admissible Fourier coefficients is combinatorially constrained by 
\ [l + (1 + %Nq) 1 I 2 ~\ < Nj < Nq + 1. In practice, that is, when the record length N is 
finite, and in particular when the signal-to-noise ratio (SNR) of the data g is small, we 
cannot demand that all the compatibility constraints be satisfied. However, checking the 
number of compatibility constraints can provide us with a confidence test on the reliability 
of the approximation fx(x). 

It is worth noting that, although we used the same notation, the set X of Eq. ()30jl can be 
different from the set X of Eq. (|18|). The former is actually the result of an algorithm acting 
on a given set of data and can be thought of as a numerical realization of the theoretical set 
X of ()18j) . A similar role is played by the numerical approximation fx(x) with respect to the 
theoretical approximation B g of Eq. PU|) . 

III. NUMERICAL EXAMPLES 

In Section [H] we illustrated a statistical procedure that allows us to split the Fourier 
coefficients into the two classes X and M and, accordingly, to write the approximation fi(x) 
[see formula (|3ip]. which, for the sake of convenience, can be rewritten as 

mi(e) _ 
m=l 

where mj(e) represents the cardinality of the set X and the eigenvalues {\k}kei (and the cor- 
responding eigenfunctions {il>k}k£i) have been suitably relabelled in a monotonic decreasing 
sequence. 
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In general, the approximation fi(x) does not satisfy the positivity constraint. How to 
incorporate effectively the positivity constraint into a regularizing scheme remains an open 
question, which has been extensively discussed in the literature (see, for instance, Refs. 
1 UM, LLii 120 and the references therein) . 
If fi(x) is not already a nonnegative function (see Fig. Q), we can aim at constructing a 
new positive approximation starting from the function fi(x) itself. The eventual negative 
part of fj(x) is mainly due to two factors: (i) the perturbation due to the noise that affects 
the coefficients {g m }mei', (h) the error due to the truncation of the series expansion. In the 
absence of additional prior information neither of these sources of error can be removed, 
but we can nevertheless look for another approximation that is nonnegative and similar to 
fi{x), according to suitable criteria. With this in mind, our task is now to seek a positive 
regularizing solution of the following type: 

f?\x) = h{x) + ^ d m ^ m {x\ (33) 

m=mx(e)+l 

where r% is an integer parameter determining the maximum number of eigenf unctions ip m {x) 
that can be used to achieve the positive solution fj\x) and the coefficients d m represent 
the unknown coefficients to be determined by requiring the function fj\x) to be nonnega- 
tive. The function fj(x) in Eq. (J33|) differs from the approximation fj(x), resulting from 
the analysis of the autocorrelation function, in the corrective finite linear combination of 
eigenfunctions ip m (x), whose purpose is to approximate and somehow compensate the errors 
leading to the negative part of fi(x). 

This kind of approach quite naturally leads us to formulate the problem of finding the 
coefficients d m as a mathematical programming problem j^jj, i.e., choosing values of a set of 
variables subject to various kinds of constraints placed on them. In particular, in our case 
we adopt a quadratic programming scheme that can be summarized as follows: Minimize 

F(d) = \\f?{x)-h{x)\\l, { _ ljiy (34) 

d = (d mi+ i, rf mi+ 2, . . . , drn x ), (35) 

subject to the constraints 

fi P \ Xl )^0, i = l,2,...,N p , (36) 



12 



where {x^Zi ls a given set of points distributed on the interval [—1, 1], and 

|A m d m - g m \ ^ e m , m = m x + 1, . . . , m x , (37) 

where e m represents an upper bound on the m th Fourier component of the noise n(y), 
computed with respect to the basis if) m . In practice only a global bound e on the noise is 
usually available, whereas the bounds e m on the single Fourier components are unknown. 
This leads us to substitute, a bit arbitrarily, e for the e m 's, thus allowing the coefficients d m 
to vary on a wider range. However, it will be shown in the numerical examples that this 
question is not a major issue, since the resulting products \ m d m generally differ from the 
corresponding g m much less than e. 

Objective function ()34j) is of a type that reflects our strategy: We look for a positive 
solution that is closest, in the sense of the L 2 norm, to the approximation fx(x), which 
is not necessarily a positive solution to our problem. Constraints (j3T?j) simply express the 
explicit requirement of positivity on a selected set of points. Regarding the constraints (|37j). 
they basically require the coefficients d m to be compatible with the Fourier coefficients of 
the data g m within the noise level. Concerning the actual numerical implementation of the 
algorithm, eigenvalues and eigenfunctions of operator A [see Eq. were computed nu- 
merically for different values of the parameter c by using the Gauss-Legendre quadrature 



22j and subsequently by diagonalizing the discretized problem by means of standard rou- 



tines [22j. The constrained optimization procedure was built around the routine E04NCF 
from the Nag Library. For every test function f(x), the corresponding data function g(y) 
was computed with Eq. ([!]) . and then noise was added. For simplicity we used only data 
corrupted by white noise simulated by computer generated random numbers uniformly dis- 
tributed in the interval [— e, e] . However, provided the assumption of independence between 
£ and £ (see Subsection III A|) . more general cases involving colored noise could be treated 



by using suitable methods such as prewhitening transformations [23j, whose discussion is 
beyond the scope of this paper. Finally, the performance of the algorithm is evaluated by 
direct comparison between the reconstructed solution and the true solution f(x). 

In Fig. ^the main steps of the analysis outlined in Subsection lH Bl are shown. The sample 
function is fi(x) = (1— x) sin 2 [4(l+a;)] with noise boundary e = 10~ 3 and corresponding SNR 
~ 40 dB. The first noiseless Fourier coefficients g^ of the image function, computed by using 
the kernel K(x, y) with c = 20, are plotted in Fig. ^\ (see the legend for numerical details). 
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Figure^3 shows the behaviour of the autocorrelation function Sg(n) along with the lines that 
indicate the statistical confidence limits we used to discriminate whether the autocorrelation 
function is substantially null. Note that the autocorrelations at n — 15, 16 have been 
correctly rejected, in spite of their quite large value, since they were abnormally inflated by 
the autocorrelations with n < 15. The approximation fx (dashed curve), computed from the 
set X obtained by analysis of the autocorrelation shown in Fig. ^3, and the actual object 
function f±(x) (solid curve) are displayed in Fig. [TJC. The excellent approximation supplied 
by fi{x) is clearly evident. Moreover, the approximation does not require any additional 
processing, since it is also nonnegative. 

For the sake of comparison, the approximation fo(x), defined by Eq. (JSJ), is reported 
in Fig. [Tp. The criterion ^ e/M with M = || /"i (a;) || z. 2 (— 1,1) gives ko(e) = 29 as trun- 
cation index. However, fo(x) computed with this value of k (not plotted) yields a very 
unsatisfactory approximation of the actual solution, and this failure can be ascribed to the 
convergence, only of weak type, of the approximation fo(x) [see Eq. (JBjl ], In this case, the 
sole constraint on the norm of the solution is not sufficient for regularizing the problem; 
therefore additional a priori information, for instance, on the first derivative of the solution, 
would be needed to achieve an acceptable approximation. Further prior information would 
lead to a smaller value of k (e) and, consequently, to good solutions, as shown in Fig. [Tp, 
where an excellent reconstruction obtained with k (e) = 27 (dashed curve) is shown. In 
contrast, with fco(e) = 28 (dashed-dotted curve) wild oscillations start appearing. 

In general, fx(x) does not satisfy the positivity constraint, in particular when the SNR 
becomes small, and so the approximation fj\x) must be computed. The basic points of 
the algorithm for constructing the positive regularized solution are summarized in Fig. 
|2]for the sample function f2(x) = exp[— (x — Xq) 2 /2a 2 ] + exp[— (x + Xq) 2 /2a 2 ] with x = 0.5 
and a = 0.1. In this example, e = 3 x 10~ 2 , and SNR ~ 6.2 dB. Since the eigenvalues of the 
operator A tend to decrease almost linearly with respect to the order index when k < <AcJ_n 
(i.e., k < 26 in this example), whereas for k > Ac/ V they go to zero exponentially fast 24], 
we expect to recover quite accurately the object-function f2{x), which is characterized as 
having the bulk of its information localized in the first values of k, even in the presence of 
such a small SNR. The analysis of the autocorrelation function 5g(n), shown in Fig. |2P, leads 
to selection of the Fourier components I = {1,3,5,7,9} with a high degree of confidence, 



since every element of X satisfies all the compatibility constraints (see Subsection IIIBj) . 
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Moreover, note that only odd components were picked up, yielding an approximation with 
the same parity of the original object function. Although the reconstruction fj, illustrated 
in Fig. EP, is quite accurate, it presents regions where it takes nonsense negative values. 
Figure EP shows the result of the constrained optimization procedure we used to achieve 
a positive solution. In this case only rnj — mj = 8 corrective components have been used, 
while the positivity constraint has been imposed on N p = 64 equidistant points of the 
interval [—1, 1]. It is worth noting that the parameter mj, which determines the number of 
corrective terms in the solution (|33J). is arbitrary and must be set manually, with the only 
condition being that the constrained optimization problem has a nonempty feasible region. 
Therefore different values of % can lead to different positive problem solutions. In this 
example we have shown the solution corresponding to the smallest value of % which gives 
a minimum of the objective function (|34|) compatible with constraints (|3l)j) and (|37j). 

Figures 01 and |U sketch the results of two further examples. In Fig. the unconstrained 
reconstruction fj(x), which represents the starting point of the optimization procedure, is 
less accurate than in the previous example. By imposing the positivity constraint we can 
see that a higher accuracy in the reconstruction of the true solution /^(x) was achieved in 
the regions where fj{x) was already positive (see the two leftmost peaks in Fig. EP). 

Finally, Fig. 0] illustrates the restoration of an edge-type object. Evidently, reconstruct- 
ing discontinuous functions within a regularizing framework based on a truncated Fourier 
expansion is a difficult task because of the Gibbs-like phenomenon. However, the overall 
regularizing procedure, including the positivity constraint on the solution, can provide an 
acceptable solution, as shown in Fig. In this case the role of the parameter % is quite 
critical. The price to pay for reaching a reasonable solution is that many corrective terms 
must be used (in this example mj = 70) and also that different feasible values of 7% can 
yield quite different reconstructions. 
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Appendix A 

Let us consider the following set of functions {fn( x )}%Lv 

/„(*) = fc*^y /p (l - \ x \r , x€[-i,i], (p^i). (A.i) 

They are nonnegative functions if x G [—1, 1] and we also have 

sup \fn(x)\=( 1 ^±±) 1/P — oo, (p^l), (A.2) 



xe[-i,i] 



2 



ll/n|| LP = l, (p^l), (A.3) 
|4f n | < const. P (A.4) 

From relation (|A.4jl it follows lim n ^oo ||A/ n || LP = 0, if p > 1. Then, for any nonnegative 
function / G L 1 (— 1, 1), we have 

IIAflLi > l\\f\\ L i , 1= mf f K(x,y)dy>0, (A.5) 

and from relation ()A.5|) it follows that the operator A~ l is continuous in the topology of the 
L l norm. 

We have thus proved that the positivity constraint restores continuity in the L 1 topology 
but not in the LP topology if p > 1 . 
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FIG. 1: Example 1: = (1 - x)sin 2 [4(l + x)\, e 



10" 



c = 20. The global SNR, defined 



as the ratio of the mean power of the noiseless data to the noise variance, was SNR ~ 40 dB. A, 
Noiseless Fourier coefficients g^. B, Modulus of the autocorrelation function. From the analysis of 
Sg(n) we have no = 14, Q = {1, 2, . . . , 13, 14}. X = {1, 2, 3, 4, 6, 7, ... , 14, 15}, each element with 
maximum inner compatibility, i.e., 13. Horizontal straight line, 95% confidence limit 1.96as(n;0) 
for a purely random sequence. This limit was used to select the elements of Q for n ^ no = 14 
(solid part). Curve, confidence limit 1.96as(n; 14) for n > no (solid part) that we used for rejecting 
the autocorrelations that are spuriously inflated by statistical fluctuations, whereas for n ^ hq 
(dashed part) it shows only how the final confidence limit 1.96as(n; 14) was reached during the 
maximization procedure for setting no [see text and, in particular, Eq. (|26|)]. C, Regularized 
solution. Solid curve, actual solution fi(x). Dashed curve, reconstruction fx(x). D, Two examples 
of the approximation fo(x) [see Eq. (JSJ)]. ko(e) was chosen as the largest integer such that ^ e/M, 
where M = \\fi(x)\\L 2 (-i,i)i that is, in the current case, k$ = 29. Solid curve, actual solution fi(x). 
Dashed curve, fo(x) computed with ko = 27 (see in particular the rightmost peak); dotted-dashed 
curve, fo(x) computed with ko = 28. The approximation fo(x) computed with ko = 29, as 
prescribed by the above truncation criterion, is not displayed since it is extremely different from 
the real solution. 
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FIG. 2: Example 2: f 2 {x) = exp[-(x - x ) 2 /2a 2 } + exp[-(x + x ) 2 /2a 2 ] with x = 0.5, a = 0.1, 
e = 3 x 10~ 2 , SNR ~ 6.2 dB, c = 20. A, Noiseless Fourier coefficients g k . B, Modulus of the 
autocorrelation function, no = 8, Q = {2, 4, 6, 8}, X = {1, 3, 5, 7, 9}. C, Unconstrained regularized 
solution. Solid curve, actual solution f2(x). Dashed curve, reconstruction fx(x). D, Comparison 

between the actual solution fa{x) and the constrained regularized solution fj\x) [see Eq. (|33j)]. 
The number of corrective eigenfunctions used was m% — mj = 8. The positivity constraint was set 
over N p = 64 points of the interval [—1,1] [see relation (|36|) ]. After the optimization procedure the 
coefficients were such that max m2 .^ m ^;m:i \d m ~ -^mdm 

I ~ 0.01 < e. 
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FIG. 3: Example 3: f 3 (x) = sin 2 [5(1 - x)} sin 2 [5(1 + x)\ [exp(-3x) + exp(-x)]. e = 5 x 10~ 2 , SNR 
~ 16.1dB, c = 20. A, Noiseless Fourier coefficients g^. B, Modulus of the autocorrelation function. 
n = 12, Q = {1,2,4,5,6,7,8,10,11,12}, Z= {1,2,3, 7,8,9,13}. Each element of Z has maximum 
inner compatibility, i.e., 6, with respect to the set Q. C, Unconstrained regularized solution. Solid 
curve, actual solution fz{x). Dashed curve, reconstruction fi{x). D, Comparison between the 

actual solution fs(x) and the constrained regularized solution f^f\x) [see Eq. (|33|) ]. The number 
of corrective eigenfunctions used was fnj — mj = 40. The positivity constraint was explicitly 
imposed on N p = 64 points of the interval [—1,1] [see relation (|36|) ], After the optimization 
procedure the coefficients were such that max^^^x \ d m — X m 9m 

I ~ 0.02 < e. 
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FIG. 4: Example 4: / 4 (x) = 1 for -0.6 < x sC 0.6 and null elsewhere, e = 8 x 10~ 2 , SNR 
~ 3.2dB, c = 20. A, Noiseless Fourier coefficients g^. B, Modulus of the autocorrelation function, 
no = 17, Q = {2, 4, 6, 17}. The set Q led to selection of the Fourier components k = 1,3, 5, 7, 18. 
However, the component at k = 18 has the minimum compatibility index, i.e., 1, which strongly 
indicates that the correlation at the lag n = 17 was spuriously generated by the noise. Then, 
X= {1,3,5,7}. C, Unconstrained regularized solution. Solid curve, actual solution fi(x). Dashed 
curve, reconstruction fj(x). D, Comparison between the actual solution fn{x) and the constrained 
regularized solution f^\ x ) [ see Eq. (|33|)] . The number of corrective eigenfunctions used was mj — 
mx = 70. The positivity constraint was set over N p = 64 points of the interval [—1,1] [see relation 
(|36j)]. After the optimization procedure the coefficients were such that inax mi ^ m ^j j d m — 
A m ffm| ^ 0.01 < e. 



